# Replication code for "Fair Share" 7-15-2020
# Ryan Brutger and Brian Rathbun
# Last modified: July 15, 2020

# This R file contains the code necessary to replicate the analysis in the main text 
# and the supplementary appendix. All of the following analyses were carried out using 
# R version 3.6.1 on a 2.8 GHz Intel Core i7 MacBookPro running macOS Mojave version 10.14.6. 

################# 1. Load libraries, load and clean data ##############

#First, set your working directory using the setwd() command:
setwd()

data <- get(load("FairShareReplicationData.RData"))

#Load required R packages (download first if necessary)
library(foreign)
library(Hmisc)
library(ggplot2)
library(psych)
library(stargazer)
library(mediation)

#Rescaling function for more easily interpretable results
rescale <- function(x){
  return((x-min(x,na.rm=TRUE))/(max(x-min(x,na.rm=TRUE),na.rm=TRUE)))
}

# Scale cooperative internationalism and  national attachment measures for easy interpretability
data$ci <- rescale(data$ci)
data$NatAtt <- rescale(data$NatAtt)

# Set critical value used for confidence intervals
c.value <-qnorm(.95)

##
## Analyze the Trade Concessions Study
## "V1" is used to indicate the data is from the Trade Concessions study

# Create data subset to those who answered the support dependent variable
dataV1 <- subset(data, is.na(data$V1_supp)==FALSE ) 
# Create data subset those who answered fairness question
dataV1_fair <- subset(dataV1[is.na(dataV1$V1_fair)==FALSE, ])

# Create Figure 1
# Calculate the mean and confidence intervals for support for each treatment

V1_T3030 <- mean(dataV1$V1_supp[dataV1$T3030==1])
se.V1_T3030 <- sqrt(var(dataV1$V1_supp[dataV1$T3030==1])/length(dataV1$V1_supp[dataV1$T3030==1])) 
ci.up.V1_T3030 <- V1_T3030 + se.V1_T3030 * c.value   
ci.low.V1_T3030 <- V1_T3030  - se.V1_T3030 * c.value 

V1_T3060 <- mean(dataV1$V1_supp[dataV1$T3060==1])
se.V1_T3060 <- sqrt(var(dataV1$V1_supp[dataV1$T3060==1])/length(dataV1$V1_supp[dataV1$T3060==1])) 
ci.up.V1_T3060 <- V1_T3060 + se.V1_T3060 * c.value   
ci.low.V1_T3060 <- V1_T3060  - se.V1_T3060 * c.value 

V1_T3090 <- mean(dataV1$V1_supp[dataV1$T3090==1])
se.V1_T3090 <- sqrt(var(dataV1$V1_supp[dataV1$T3090==1])/length(dataV1$V1_supp[dataV1$T3090==1])) 
ci.up.V1_T3090 <- V1_T3090 + se.V1_T3090 * c.value   
ci.low.V1_T3090 <- V1_T3090  - se.V1_T3090 * c.value

V1_T6030 <- mean(dataV1$V1_supp[dataV1$T6030==1])
se.V1_T6030 <- sqrt(var(dataV1$V1_supp[dataV1$T6030==1])/length(dataV1$V1_supp[dataV1$T6030==1])) 
ci.up.V1_T6030 <- V1_T6030 + se.V1_T6030 * c.value   
ci.low.V1_T6030 <- V1_T6030  - se.V1_T6030 * c.value

V1_T9030 <- mean(dataV1$V1_supp[dataV1$T9030==1])
se.V1_T9030 <- sqrt(var(dataV1$V1_supp[dataV1$T9030==1])/length(dataV1$V1_supp[dataV1$T9030==1])) 
ci.up.V1_T9030 <- V1_T9030 + se.V1_T9030 * c.value   
ci.low.V1_T9030 <- V1_T9030  - se.V1_T9030 * c.value

V1_T6060 <- mean(dataV1$V1_supp[dataV1$T6060==1])
se.V1_T6060 <- sqrt(var(dataV1$V1_supp[dataV1$T6060==1])/length(dataV1$V1_supp[dataV1$T6060==1])) 
ci.up.V1_T6060 <- V1_T6060 + se.V1_T6060 * c.value   
ci.low.V1_T6060 <- V1_T6060  - se.V1_T6060 * c.value

V1_T6090 <- mean(dataV1$V1_supp[dataV1$T6090==1])
se.V1_T6090 <- sqrt(var(dataV1$V1_supp[dataV1$T6090==1])/length(dataV1$V1_supp[dataV1$T6090==1])) 
ci.up.V1_T6090 <- V1_T6090 + se.V1_T6090 * c.value   
ci.low.V1_T6090 <- V1_T6090  - se.V1_T6090 * c.value

V1_T9090 <- mean(dataV1$V1_supp[dataV1$T9090==1])
se.V1_T9090 <- sqrt(var(dataV1$V1_supp[dataV1$T9090==1])/length(dataV1$V1_supp[dataV1$T9090==1])) 
ci.up.V1_T9090 <- V1_T9090 + se.V1_T9090 * c.value   
ci.low.V1_T9090 <- V1_T9090  - se.V1_T9090 * c.value

V1_T9060 <- mean(dataV1$V1_supp[dataV1$T9060==1])
se.V1_T9060 <- sqrt(var(dataV1$V1_supp[dataV1$T9060==1])/length(dataV1$V1_supp[dataV1$T9060==1])) 
ci.up.V1_T9060 <- V1_T9060 + se.V1_T9060 * c.value   
ci.low.V1_T9060 <- V1_T9060  - se.V1_T9060 * c.value

# Consildate means and confidence intervals for plot 
V1.means.abs <- c(V1_T3030, V1_T3060, V1_T3090, V1_T6030, V1_T6060, V1_T6090, V1_T9030, V1_T9060, V1_T9090)
V1.ci.up.abs <- c(ci.up.V1_T3030, ci.up.V1_T3060, ci.up.V1_T3090, ci.up.V1_T6030, ci.up.V1_T6060, ci.up.V1_T6090, ci.up.V1_T9030, ci.up.V1_T9060, ci.up.V1_T9090)
V1.ci.low.abs <- c(ci.low.V1_T3030, ci.low.V1_T3060, ci.low.V1_T3090, ci.low.V1_T6030, ci.low.V1_T6060, ci.low.V1_T6090, ci.low.V1_T9030, ci.low.V1_T9060, ci.low.V1_T9090)
V1.names.abs <- c("30-30","30-60", "30-90", "60-30", "60-60","60-90", "90-30", "90-60", "90-90")

## Generate plot from smallest to largest tariff reductions (adjust width/height as necessary)
par(mar=c(3,5,1,2))
plot(1:9, V1.means.abs , main = "", 
     ylab = "", xlab = "", cex.lab=2,
     ylim = c(-.35, .9), xlim = c(0.75, 9.25), pch = 19, type = "p", xaxt="n", cex=1)
mtext("Average Support Score", side=2, at=.25, cex=1.9, line=2.5)
abline(h=0, col="blue", lty=2)
for(i in 1:length(V1.means.abs)){
  lines(c(i,i), c(V1.ci.up.abs[i], V1.ci.low.abs[i]))
  mtext(V1.names.abs[i], side=1, at=i, 1.12, cex=1)
}

# Save the preceeding estimates by equal, favorable, and unfavorable treatments
# which will be used to generate Figure 3

# The following 3 lines are for the equal treatment plot
V1.means.eq <- c(V1_T3030, V1_T6060, V1_T9090)
V1.ci.up.eq <- c(ci.up.V1_T3030, ci.up.V1_T6060, ci.up.V1_T9090)
V1.ci.low.eq <- c(ci.low.V1_T3030, ci.low.V1_T6060, ci.low.V1_T9090)

# The following 3 lines are for the unfavorable treatment plot
V1.means.un <- c(V1_T9030, V1_T9060, V1_T6030)
V1.ci.up.un <- c(ci.up.V1_T9030, ci.up.V1_T9060, ci.up.V1_T6030)
V1.ci.low.un <- c(ci.low.V1_T9030, ci.low.V1_T9060, ci.low.V1_T6030)

# The following 3 lines are for the favorable treatment plot
V1.means.fav <- c(V1_T3090, V1_T3060, V1_T6090)
V1.ci.up.fav <- c(ci.up.V1_T3090, ci.up.V1_T3060, ci.up.V1_T6090)
V1.ci.low.fav <- c(ci.low.V1_T3090, ci.low.V1_T3060, ci.low.V1_T6090)

# Now create Figure 2
# Analyze Trade Concession treatments by equal, unfavorable, or favorable tariff reductions 

# Begin with the means and confidence intervals for support
V1_equal <- mean(dataV1$V1_supp[dataV1$V1_equal==1])
se.V1_equal <- sqrt(var(dataV1$V1_supp[dataV1$V1_equal==1])/length(dataV1$V1_supp[dataV1$V1_equal==1])) 
ci.up.V1_equal <- V1_equal + se.V1_equal * c.value   
ci.low.V1_equal <- V1_equal  - se.V1_equal * c.value 

V1_unfav <- mean(dataV1$V1_supp[dataV1$V1_unfav==1])
se.V1_unfav <- sqrt(var(dataV1$V1_supp[dataV1$V1_unfav==1])/length(dataV1$V1_supp[dataV1$V1_unfav==1])) 
ci.up.V1_unfav <- V1_unfav + se.V1_unfav * c.value   
ci.low.V1_unfav <- V1_unfav  - se.V1_unfav * c.value

V1_fav <- mean(dataV1$V1_supp[dataV1$V1_fav==1])
se.V1_fav <- sqrt(var(dataV1$V1_supp[dataV1$V1_fav==1])/length(dataV1$V1_supp[dataV1$V1_fav==1])) 
ci.up.V1_fav <- V1_fav + se.V1_fav * c.value   
ci.low.V1_fav <- V1_fav  - se.V1_fav * c.value

# Consolidate support estimates for plot
V1.means.s <- c(V1_equal, V1_unfav, V1_fav)
V1.ci.up.s <- c(ci.up.V1_equal, ci.up.V1_unfav, ci.up.V1_fav)
V1.ci.low.s <- c(ci.low.V1_equal, ci.low.V1_unfav, ci.low.V1_fav)
V1.names.s <- c("Equal", "Unfavorable", "Favorable")

## Now do the same for fairness

V1_equal <- mean(dataV1_fair$V1_fair[dataV1_fair$V1_equal==1])
se.V1_equal <- sqrt(var(dataV1_fair$V1_fair[dataV1_fair$V1_equal==1])/length(dataV1_fair$V1_fair[dataV1_fair$V1_equal==1])) 
ci.up.V1_equal <- V1_equal + se.V1_equal * c.value   
ci.low.V1_equal <- V1_equal  - se.V1_equal * c.value 

V1_unfav <- mean(dataV1_fair$V1_fair[dataV1_fair$V1_unfav==1])
se.V1_unfav <- sqrt(var(dataV1_fair$V1_fair[dataV1_fair$V1_unfav==1])/length(dataV1_fair$V1_fair[dataV1_fair$V1_unfav==1])) 
ci.up.V1_unfav <- V1_unfav + se.V1_unfav * c.value   
ci.low.V1_unfav <- V1_unfav  - se.V1_unfav * c.value

V1_fav <- mean(dataV1_fair$V1_fair[dataV1_fair$V1_fav==1])
se.V1_fav <- sqrt(var(dataV1_fair$V1_fair[dataV1_fair$V1_fav==1])/length(dataV1_fair$V1_fair[dataV1_fair$V1_fav==1])) 
ci.up.V1_fav <- V1_fav + se.V1_fav * c.value   
ci.low.V1_fav <- V1_fav  - se.V1_fav * c.value

# Consolidate fairness estimates for plot
V1.means.f <- c(V1_equal, V1_unfav, V1_fav)
V1.ci.up.f <- c(ci.up.V1_equal, ci.up.V1_unfav, ci.up.V1_fav)
V1.ci.low.f <- c(ci.low.V1_equal, ci.low.V1_unfav, ci.low.V1_fav)
V1.names.f <- c("Equal", "Unfavorable", "Favorable")

# Generate Figure 2 (adjust width/height as necessary)
par(mar=c(3,6,1,2))
plot(1:3, V1.means.f , main = "", 
     ylab = "", xlab = "", cex.lab=2,
     ylim = c(-.35, .9), xlim = c(0.75, 3.25), pch = 19, type = "p", xaxt="n", cex=1)
mtext("Average Fairness \n & Support Scores  ", side=2, at=.25, cex=1.9, line=2.5)
abline(h=0, col="blue", lty=2)
for(i in 1:length(V1.means.f)){
  lines(c(i,i), c(V1.ci.up.f[i], V1.ci.low.f[i]))
  mtext(V1.names.f[i], side=1, at=i, 1.12, cex=1.9)
}
points(1.15:3.15, V1.means.s, pch=22, cex=1.5)
for(i in 1:length(V1.means.s)){
  lines(c(i+.15,i+.15), c(V1.ci.up.s[i], V1.ci.low.s[i]))
}
text(3.01, .85, labels = "Fairness", cex=1.5)
text(3, .75, labels = "Support", cex=1.5)
points(2.7, .85, pch=19, cex=1.5)
points(2.7, .75, pch=22, cex=1.5)

# Analysis for quantities of interest cited in the text relating to Figures 2 and 3

#Create measure of those expressing positve support (or not)
dataV1$dich_supp <- ifelse(dataV1$V1_supp>0, 1, 0)

# Generate percent supporting the equal treatment and the difference between that and the favorable treatment
eq_fav <- t.test(dataV1$dich_supp[dataV1$V1_equal==1], dataV1$dich_supp[dataV1$V1_fav==1])
eq_fav$estimate[1] # support for equal: 51%
eq_fav$estimate[1] - eq_fav$estimate[2] #-0.01583067
eq_fav$p.value #p=0.47

# Generate support score difference and  between the equal and the favorable treatment for footnote 7
eq_fav <- t.test(dataV1$V1_supp[dataV1$V1_equal==1], dataV1$V1_supp[dataV1$V1_fav==1])
eq_fav$estimate[1] - eq_fav$estimate[2] #-0.04
eq_fav$p.value #p=0.32

# Generate the difference in the percent supporting  between the equal and the unfavorable treatments
eq_unfav <- t.test(dataV1$dich_supp[dataV1$V1_equal==1], dataV1$dich_supp[dataV1$V1_unfav==1])
eq_unfav$estimate[1] - eq_unfav$estimate[2] #20%
eq_unfav$p.value #p<0.01

# Generate the difference in the percent supporting  between the favorable and the unfavorable treatments
eq_unfav <- t.test(dataV1$dich_supp[dataV1$V1_fav==1], dataV1$dich_supp[dataV1$V1_unfav==1])
eq_unfav$estimate[1] - eq_unfav$estimate[2] #21%
eq_unfav$p.value #p<0.01

# Analysis for footnote 8 and the note to Figure 3, which shows the lowest p-value 
# in the favorable condition for fairness is p=0.34
V13060_3090 <- t.test(dataV1_fair$V1_fair[dataV1_fair$T3060==1], dataV1_fair$V1_fair[dataV1_fair$T3090==1])
V13060_3090$estimate[1] - V13060_3090$estimate[2] #-0.02
V13060_3090$p.value #p<0.81

V13060_6090 <- t.test(dataV1_fair$V1_fair[dataV1_fair$T3060==1], dataV1_fair$V1_fair[dataV1_fair$T6090==1])
V13060_6090$estimate[1] - V13060_6090$estimate[2] #-0.07
V13060_6090$p.value #p=0.34

V13090_6090 <- t.test(dataV1_fair$V1_fair[dataV1_fair$T3090==1], dataV1_fair$V1_fair[dataV1_fair$T6090==1])
V13090_6090$estimate[1] - V13090_6090$estimate[2] #-0.05
V13090_6090$p.value #p<0.48

# Analysis for discussion of panel "c" of Figure 3

# Change in fairness from 60/30 to 90/30
V19030_6030 <- t.test(dataV1_fair$V1_fair[dataV1_fair$T9030==1], dataV1_fair$V1_fair[dataV1_fair$T6030==1])
V19030_6030$estimate[1] - V19030_6030$estimate[2] #-0.23
V19030_6030$p.value #p<0.01

# Change in fairness from 90/60 to 60/30
V16030_9060 <- t.test(dataV1_fair$V1_fair[dataV1_fair$T6030==1], dataV1_fair$V1_fair[dataV1_fair$T9060==1])
V16030_9060$estimate[2] - V16030_9060$estimate[1] #-0.12
V16030_9060$p.value #p=0.10

# Change in fairness from 90/60 to 90/30
V19030_9060 <- t.test(dataV1_fair$V1_fair[dataV1_fair$T9030==1], dataV1_fair$V1_fair[dataV1_fair$T9060==1])
V19030_9060$estimate[1] - V19030_9060$estimate[2] #-0.10
V19030_9060$p.value #p<0.18

# Analysis for footnote 9 and note to Figure 3, which shows none of the p-values for differences in the 
# support score amongst the favorable treatments reaches traditional levels of significance

# Change in support from 30/90 to 30/60
V13060_3090 <- t.test(dataV1$V1_supp[dataV1$T3060==1], dataV1$V1_supp[dataV1$T3090==1])
V13060_3090$estimate[1] - V13060_3090$estimate[2] #-0.07
V13060_3090$p.value #p=0.31

# Change in support from 30/60 to 60/90
V13060_6090 <- t.test(dataV1$V1_supp[dataV1$T6090==1], dataV1$V1_supp[dataV1$T3060==1])
V13060_6090$estimate[1] - V13060_6090$estimate[2]  #0.11
V13060_6090$p.value #p=0.13

# Change in support from 60/90 to 30/90
V13090_6090 <- t.test(dataV1$V1_supp[dataV1$T3090==1], dataV1$V1_supp[dataV1$T6090==1])
V13090_6090$estimate[1] - V13090_6090$estimate[2] #-0.04
V13090_6090$p.value #p=0.62


# Create Figure 3
# This requires generating three plots for panels a, b, and c of the figure
# The support measures are already saved from above, so now proceed to do the same for the fairness measure

V1_T3030 <- mean(dataV1_fair$V1_fair[dataV1_fair$T3030==1])
se.V1_T3030 <- sqrt(var(dataV1_fair$V1_fair[dataV1_fair$T3030==1])/length(dataV1_fair$V1_fair[dataV1_fair$T3030==1])) 
ci.up.V1_T3030 <- V1_T3030 + se.V1_T3030 * c.value   
ci.low.V1_T3030 <- V1_T3030  - se.V1_T3030 * c.value 

V1_T3060 <- mean(dataV1_fair$V1_fair[dataV1_fair$T3060==1])
se.V1_T3060 <- sqrt(var(dataV1_fair$V1_fair[dataV1_fair$T3060==1])/length(dataV1_fair$V1_fair[dataV1_fair$T3060==1])) 
ci.up.V1_T3060 <- V1_T3060 + se.V1_T3060 * c.value   
ci.low.V1_T3060 <- V1_T3060  - se.V1_T3060 * c.value 

V1_T3090 <- mean(dataV1_fair$V1_fair[dataV1_fair$T3090==1])
se.V1_T3090 <- sqrt(var(dataV1_fair$V1_fair[dataV1_fair$T3090==1])/length(dataV1_fair$V1_fair[dataV1_fair$T3090==1])) 
ci.up.V1_T3090 <- V1_T3090 + se.V1_T3090 * c.value   
ci.low.V1_T3090 <- V1_T3090  - se.V1_T3090 * c.value

V1_T6030 <- mean(dataV1_fair$V1_fair[dataV1_fair$T6030==1])
se.V1_T6030 <- sqrt(var(dataV1_fair$V1_fair[dataV1_fair$T6030==1])/length(dataV1_fair$V1_fair[dataV1_fair$T6030==1])) 
ci.up.V1_T6030 <- V1_T6030 + se.V1_T6030 * c.value   
ci.low.V1_T6030 <- V1_T6030  - se.V1_T6030 * c.value

V1_T9030 <- mean(dataV1_fair$V1_fair[dataV1_fair$T9030==1])
se.V1_T9030 <- sqrt(var(dataV1_fair$V1_fair[dataV1_fair$T9030==1])/length(dataV1_fair$V1_fair[dataV1_fair$T9030==1])) 
ci.up.V1_T9030 <- V1_T9030 + se.V1_T9030 * c.value   
ci.low.V1_T9030 <- V1_T9030  - se.V1_T9030 * c.value

V1_T6060 <- mean(dataV1_fair$V1_fair[dataV1_fair$T6060==1])
se.V1_T6060 <- sqrt(var(dataV1_fair$V1_fair[dataV1_fair$T6060==1])/length(dataV1_fair$V1_fair[dataV1_fair$T6060==1])) 
ci.up.V1_T6060 <- V1_T6060 + se.V1_T6060 * c.value   
ci.low.V1_T6060 <- V1_T6060  - se.V1_T6060 * c.value

V1_T6090 <- mean(dataV1_fair$V1_fair[dataV1_fair$T6090==1])
se.V1_T6090 <- sqrt(var(dataV1_fair$V1_fair[dataV1_fair$T6090==1])/length(dataV1_fair$V1_fair[dataV1_fair$T6090==1])) 
ci.up.V1_T6090 <- V1_T6090 + se.V1_T6090 * c.value   
ci.low.V1_T6090 <- V1_T6090  - se.V1_T6090 * c.value

V1_T9090 <- mean(dataV1_fair$V1_fair[dataV1_fair$T9090==1])
se.V1_T9090 <- sqrt(var(dataV1_fair$V1_fair[dataV1_fair$T9090==1])/length(dataV1_fair$V1_fair[dataV1_fair$T9090==1])) 
ci.up.V1_T9090 <- V1_T9090 + se.V1_T9090 * c.value   
ci.low.V1_T9090 <- V1_T9090  - se.V1_T9090 * c.value

V1_T9060 <- mean(dataV1_fair$V1_fair[dataV1_fair$T9060==1])
se.V1_T9060 <- sqrt(var(dataV1_fair$V1_fair[dataV1_fair$T9060==1])/length(dataV1_fair$V1_fair[dataV1_fair$T9060==1])) 
ci.up.V1_T9060 <- V1_T9060 + se.V1_T9060 * c.value   
ci.low.V1_T9060 <- V1_T9060  - se.V1_T9060 * c.value

## Save the fairness measures for the equal treatment
V1.means.f.eq <- c(V1_T3030, V1_T6060, V1_T9090)
V1.ci.up.f.eq <- c(ci.up.V1_T3030, ci.up.V1_T6060, ci.up.V1_T9090)
V1.ci.low.f.eq <- c(ci.low.V1_T3030, ci.low.V1_T6060, ci.low.V1_T9090)
V1.names.f.eq <- c("30-30","60-60", "90-90")

#  Save the fairness measures for the unfavorable treatment 
V1.means.f.un <- c(V1_T9030, V1_T9060, V1_T6030)
V1.ci.up.f.un <- c(ci.up.V1_T9030, ci.up.V1_T9060, ci.up.V1_T6030)
V1.ci.low.f.un <- c(ci.low.V1_T9030, ci.low.V1_T9060, ci.low.V1_T6030)
V1.names.f.un <- c("90-30","90-60", "60-30")

# Save the fairness measures for the favorable treatment plot
V1.means.f.fav <- c(V1_T3090, V1_T3060, V1_T6090)
V1.ci.up.f.fav <- c(ci.up.V1_T3090, ci.up.V1_T3060, ci.up.V1_T6090)
V1.ci.low.f.fav <- c(ci.low.V1_T3090, ci.low.V1_T3060, ci.low.V1_T6090)
V1.names.f.fav <- c("30-90","30-60", "60-90")

# Generate panel "a" of Figure 3
par(mar=c(3,5,1,2))
plot(1:3, V1.means.f.eq , main = "", 
     ylab = "", xlab = "", cex.lab=2,
     ylim = c(-.05, .85), xlim = c(0.75, 3.25), pch = 19, type = "p", xaxt="n", cex=1)
mtext("Average Fairness & Support Scores", side=2, at=.35, cex=1.8, line=2.5)
abline(h=0, col="blue", lty=2)
for(i in 1:length(V1.means.f)){
  lines(c(i,i), c(V1.ci.up.f.eq[i], V1.ci.low.f.eq[i]))
  mtext(V1.names.f.eq[i], side=1, at=i, 1.12, cex=1)
}
points(1.15:3.15, V1.means.eq, pch=22, cex=1.5)
for(i in 1:length(V1.means.eq)){
  lines(c(i+.15,i+.15), c(V1.ci.up.eq[i], V1.ci.low.eq[i]))
}
text(3, .8, labels = "Fairness", cex=1.5)
text(3, .73, labels = "Support", cex=1.5)
points(2.7, .8, pch=19, cex=1.5)
points(2.7, .73, pch=22, cex=1.5)

# Generate panel "b" of Figure 3
par(mar=c(3,5,1,2))
plot(1:3, V1.means.f.fav , main = "", 
     ylab = "", xlab = "", cex.lab=2,
     ylim = c(-.35, .85), xlim = c(0.75, 3.25), pch = 19, type = "p", xaxt="n", cex=1)
mtext("", side=2, at=.4, cex=1.8, line=2.5)
abline(h=0, col="blue", lty=2)
for(i in 1:length(V1.means.f.fav)){
  lines(c(i,i), c(V1.ci.up.f.fav[i], V1.ci.low.f.fav[i]))
  mtext(V1.names.f.eq[i], side=1, at=i, 1.12, cex=1)
}
points(1.15:3.15, V1.means.fav, pch=22, cex=1.5)
for(i in 1:length(V1.means.fav)){
  lines(c(i+.15,i+.15), c(V1.ci.up.fav[i], V1.ci.low.fav[i]))
}
text(3, .84, labels = "Fairness", cex=1.5)
text(3, .75, labels = "Support", cex=1.5)
points(2.7, .84, pch=19, cex=1.5)
points(2.7, .75, pch=22, cex=1.5)

# Generate panel "c" of Figure 3
par(mar=c(3,5,1,2))
plot(1:3, V1.means.f.un , main = "", 
     ylab = "", xlab = "", cex.lab=2,
     ylim = c(-.35, .85), xlim = c(0.75, 3.25), pch = 19, type = "p", xaxt="n", cex=1)
mtext("Average Fairness & Support Scores", side=2, at=.2, cex=1.8, line=2.5)
abline(h=0, col="blue", lty=2)
for(i in 1:length(V1.means.f.un)){
  lines(c(i,i), c(V1.ci.up.f.un[i], V1.ci.low.f.un[i]))
  mtext(V1.names.f.eq[i], side=1, at=i, 1.12, cex=1)
}
points(1.15:3.15, V1.means.un, pch=22, cex=1.5)
for(i in 1:length(V1.means.un)){
  lines(c(i+.15,i+.15), c(V1.ci.up.un[i], V1.ci.low.un[i]))
}
text(3, .84, labels = "Fairness", cex=1.5)
text(3, .75, labels = "Support", cex=1.5)
points(2.7, .84, pch=19, cex=1.5)
points(2.7, .75, pch=22, cex=1.5)

# Create Figure 4
# Figure uses the seperate figures for panels a, b, and c

# Calculations for panel "a" of Figure 4: mediation analysis comparing equal versus favorable agreements

# subset to equal and favorable treatment
set.seed(43216)
dataV1eq.fav <- subset(dataV1[dataV1$V1_equal==1 | dataV1$V1_fav==1, ])

## fairness of agreement as a mediator
fair2.med <- lm(V1_fair ~ V1_equal + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1eq.fav)
fair2.dv <- lm(V1_supp ~ V1_fair + V1_equal + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1eq.fav)
fair2.med2 <- mediate(fair2.med, fair2.dv, treat="V1_equal", mediator="V1_fair", sims=1500)

# Extract quantities of interest for plot (note: plots can made using the plot function, but they are displayed horizontally by default)
fair2.effects <- c(fair2.med2[1], fair2.med2[9], fair2.med2$tau.coef)
fair2.ci.low <- c(fair2.med2$d1.ci[1], fair2.med2$z0.ci[1], fair2.med2$tau.ci[1])
fair2.ci.high <- c(fair2.med2$d1.ci[2], fair2.med2$z0.ci[2], fair2.med2$tau.ci[2])

fair2.med.names <- c("ACME \n", "ADE \n ", "Total \n Effect")

## Generate for panel "a" of Figure 4
par(mar=c(3,5,1,2))
plot(1:3, fair2.effects , main = "", 
     ylab = "", xlab = "", cex.lab=2,
     ylim = c(-.35, .65), xlim = c(0.75, 3.25), pch = 19, type = "p", xaxt="n", cex=1)
mtext("", side=2, at=.2, cex=1.7, line=2.5)
abline(h=0, col="blue", lty=2)
for(i in 1:length(fair2.med.names)){
  lines(c(i,i), c(fair2.ci.high[i], fair2.ci.low[i]))
  mtext(fair2.med.names[i], side=1, at=i, 1.75, cex=1.3)
}

# Calculations for panel "b" of Figure 4: mediation analysis comparing equal versus unfavorable agreements

# subset to equal and unfavorable treatment
set.seed(43216)
dataV1eq.unfav <- subset(dataV1[dataV1$V1_equal==1 | dataV1$V1_unfav==1, ])

## fairness of agreement as a mediator
fair.med <- lm(V1_fair ~ V1_equal + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1eq.unfav)
fair.dv <- lm(V1_supp ~ V1_fair + V1_equal + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1eq.unfav)
fair.med2 <- mediate(fair.med, fair.dv, treat="V1_equal", mediator="V1_fair", sims=1500)

# Extract quantities of interest for plot (note: plots can made using the plot function, but they are displayed horizontally by default)
fair.effects <- c(fair.med2[1], fair.med2[9], fair.med2$tau.coef)
fair.ci.low <- c(fair.med2$d1.ci[1], fair.med2$z0.ci[1], fair.med2$tau.ci[1])
fair.ci.high <- c(fair.med2$d1.ci[2], fair.med2$z0.ci[2], fair.med2$tau.ci[2])

fair.med.names <- c("ACME \n", "ADE \n ", "Total \n Effect")

## Generate for panel "b" of Figure 4
par(mar=c(3,5,1,2))
plot(1:3, fair.effects , main = "", 
     ylab = "", xlab = "", cex.lab=2,
     ylim = c(-.35, .65), xlim = c(0.75, 3.25), pch = 19, type = "p", xaxt="n", cex=1)
mtext("", side=2, at=.2, cex=1.7, line=2.5)
abline(h=0, col="blue", lty=2)
for(i in 1:length(fair.med.names)){
  lines(c(i,i), c(fair.ci.high[i], fair.ci.low[i]))
  mtext(fair.med.names[i], side=1, at=i, 1.75, cex=1.3)
}

# For quantities of interest cited in the text
summary(fair.med2) 
# ACME 0.51
# p < 0.01  

# Calculations for panel "c" of Figure 4: mediation analysis comparing favorable versus unfavorable agreements

# subset to favorable and unfavorable treatment
set.seed(43216)
dataV1fav.unfav <- subset(dataV1[dataV1$V1_unfav==1 | dataV1$V1_fav==1, ])

## fairness of agreement as a mediator
fair3.med <- lm(V1_fair ~ V1_fav + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1fav.unfav)
fair3.dv <- lm(V1_supp ~ V1_fair + V1_fav + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1fav.unfav)
fair3.med2 <- mediate(fair3.med, fair3.dv, treat="V1_fav", mediator="V1_fair", sims=1500)

# Extract quantities of interest for plot (note: plots can made using the plot function, but they are displayed horizontally by default)
fair3.effects <- c(fair3.med2[1], fair3.med2[9], fair3.med2$tau.coef)
fair3.ci.low <- c(fair3.med2$d1.ci[1], fair3.med2$z0.ci[1], fair3.med2$tau.ci[1])
fair3.ci.high <- c(fair3.med2$d1.ci[2], fair3.med2$z0.ci[2], fair3.med2$tau.ci[2])

fair3.med.names <- c("ACME \n", "ADE \n ", "Total \n Effect")

## Generate for panel "c" of Figure 4
par(mar=c(3,5,1,2))
plot(1:3, fair3.effects , main = "", 
     ylab = "", xlab = "", cex.lab=2,
     ylim = c(-.35, .65), xlim = c(0.75, 3.25), pch = 19, type = "p", xaxt="n", cex=1)
mtext("", side=2, at=.2, cex=1.7, line=2.5)
abline(h=0, col="blue", lty=2)
for(i in 1:length(fair3.med.names)){
  lines(c(i,i), c(fair3.ci.high[i], fair3.ci.low[i]))
  mtext(fair3.med.names[i], side=1, at=i, 1.75, cex=1.3)
}

# For quantities of interest cited in the text
summary(fair3.med2) 
# ACME 0.33          
# p <0.01
# proportion mediated 0.59 (about 60 percent)

##
## Begin analysis of Trade Balance study
##

# Create subset to those who answered support for renegotiation dependent variable
# "V3" is used to identify data from the Trade Balance study
dataV3 <- subset(data, is.na(data$V3_App_reneg)==FALSE )

# Analysis for Figure 5
# Calculate means and confidence intervals for renegotiation support for favorable, unfavorable, and balanced treatments
V3_fav <- mean(dataV3$V3_App_reneg[dataV3$V3_fav==1])
se.V3_fav <- sqrt(var(dataV3$V3_App_reneg[dataV3$V3_fav==1])/length(dataV3$V3_App_reneg[dataV3$V3_fav==1])) 
ci.up.V3_fav <- V3_fav + se.V3_fav * c.value   
ci.low.V3_fav <- V3_fav  - se.V3_fav * c.value 

V3_unfav <- mean(dataV3$V3_App_reneg[dataV3$V3_unfav==1])
se.V3_unfav <- sqrt(var(dataV3$V3_App_reneg[dataV3$V3_unfav==1])/length(dataV3$V3_App_reneg[dataV3$V3_unfav==1])) 
ci.up.V3_unfav <- V3_unfav + se.V3_unfav * c.value   
ci.low.V3_unfav <- V3_unfav  - se.V3_unfav * c.value 

V3_bal <- mean(dataV3$V3_App_reneg[dataV3$V3_bal==1])
se.V3_bal <- sqrt(var(dataV3$V3_App_reneg[dataV3$V3_bal==1])/length(dataV3$V3_App_reneg[dataV3$V3_bal==1])) 
ci.up.V3_bal <- V3_bal + se.V3_bal * c.value   
ci.low.V3_bal <- V3_bal  - se.V3_bal * c.value 

# Save consolidated estimates for plot
V3.means.r <- c(V3_bal, V3_unfav, V3_fav)
V3.ci.up.r <- c(ci.up.V3_bal, ci.up.V3_unfav, ci.up.V3_fav)
V3.ci.low.r <- c(ci.low.V3_bal, ci.low.V3_unfav, ci.low.V3_fav)
V3.names.r <- c("Balanced", "Unfavorable", "Favorable")

# Calculate means and confidence intervals for fairness for favorable, unfavorable, and balanced treatments
# subset to those who answered fairness dependent variable
dataV3 <- subset(data, is.na(data$V3_fair)==FALSE )

V3_fav <- mean(dataV3$V3_fair[dataV3$V3_fav==1])
se.V3_fav <- sqrt(var(dataV3$V3_fair[dataV3$V3_fav==1])/length(dataV3$V3_fair[dataV3$V3_fav==1])) 
ci.up.V3_fav <- V3_fav + se.V3_fav * c.value   
ci.low.V3_fav <- V3_fav  - se.V3_fav * c.value 

V3_unfav <- mean(dataV3$V3_fair[dataV3$V3_unfav==1])
se.V3_unfav <- sqrt(var(dataV3$V3_fair[dataV3$V3_unfav==1])/length(dataV3$V3_fair[dataV3$V3_unfav==1])) 
ci.up.V3_unfav <- V3_unfav + se.V3_unfav * c.value   
ci.low.V3_unfav <- V3_unfav  - se.V3_unfav * c.value 

V3_bal <- mean(dataV3$V3_fair[dataV3$V3_bal==1])
se.V3_bal <- sqrt(var(dataV3$V3_fair[dataV3$V3_bal==1])/length(dataV3$V3_fair[dataV3$V3_bal==1])) 
ci.up.V3_bal <- V3_bal + se.V3_bal * c.value   
ci.low.V3_bal <- V3_bal  - se.V3_bal * c.value 

V3.means <- c(V3_bal, V3_unfav, V3_fav)
V3.ci.up <- c(ci.up.V3_bal, ci.up.V3_unfav, ci.up.V3_fav)
V3.ci.low <- c(ci.low.V3_bal, ci.low.V3_unfav, ci.low.V3_fav)
V3.names <- c("Balanced", "Unfavorable", "Favorable")

# Generate Figure 5 (adjust width/height as necessary)
par(mar=c(3,6,1,2))
plot(1:3, V3.means , main = "", 
     ylab = "", xlab = "", cex.lab=2,
     ylim = c(-.25, 0.95), xlim = c(0.75, 3.25), pch = 19, type = "p", xaxt="n", cex=1.5)
mtext("Average Fairness \n and Renogitiation Scores", side=2, at=.25, cex=1.7, line=2.45)
abline(h=0, col="blue", lty=2)
for(i in 1:length(V3.means)){
  lines(c(i,i), c(V3.ci.up[i], V3.ci.low[i]))
  mtext(V3.names[i], side=1, at=i, 1.12, cex=1.7)
}
points(1.05:3.05, V3.means.r, pch=22, cex=1.5)
for(i in 1:length(V3.means)){
  lines(c(i+.05,i+.05), c(V3.ci.up.r[i], V3.ci.low.r[i]))
}
text(2.27, .86, labels = "Fairness", cex=1.5)
text(2.68, .95, labels = "Supp. for Renegotiation", cex=1.5)
points(1.98, .86, pch=19, cex=1.5)
points(1.98, .95, pch=22, cex=1.5)

# Analysis for quantities of interest cited in the text from the discussion of Figure 5

# Change in fairness from unfavorable to equal (balanced) treatment
eq_unfav <- t.test(dataV3$V3_fair[dataV3$V3_bal==1], dataV3$V3_fair[dataV3$V3_unfav==1])
eq_unfav$estimate[1] - eq_unfav$estimate[2] #0.43
eq_unfav$p.value #p<0.01

# Change in fairness from balanced to favorable treatment
fav_eq <- t.test(dataV3$V3_fair[dataV3$V3_fav==1], dataV3$V3_fair[dataV3$V3_bal==1])
fav_eq$estimate[2] - fav_eq$estimate[1] #0.14
fav_eq$p.value #p<0.01

# Change in percent who believe agreement is fair
# create dichotomous  measure
dataV3$dich_fair <- ifelse(dataV3$V3_fair>0, 1, 0)

# Change in percent who believe agreement is fair from unfavorable to balanced
eq_unfav <- t.test(dataV3$dich_fair[dataV3$V3_bal==1], dataV3$dich_fair[dataV3$V3_unfav==1])
eq_unfav$estimate[1] - eq_unfav$estimate[2] #0.13%
eq_unfav$p.value #p<0.01

# Change in percent who believe agreement is fair from favorable to balanced
fav_eq <- t.test(dataV3$dich_fair[dataV3$V3_bal==1], dataV3$dich_fair[dataV3$V3_fav==1])
fav_eq$estimate[1] - fav_eq$estimate[2]  #7%
fav_eq$p.value #p<0.01

# Now quntities of interest for support for renegotiation
dataV3 <- subset(data, is.na(data$V3_App_reneg)==FALSE )

# Change in support for rengotiation from favorable to balanced 
fav_eq <- t.test(dataV3$V3_App_reneg[dataV3$V3_bal==1], dataV3$V3_App_reneg[dataV3$V3_fav==1])
fav_eq$estimate[1] - fav_eq$estimate[2] #0.07
fav_eq$p.value #p=0.08

# create dichotomous measure for percent supporting renegotiation
dataV3$dich_reneg <- ifelse(dataV3$V3_App_reneg>0, 1, 0)

# Change in percent supporting renegotiation from balanced to favorable
fav_eq <- t.test(dataV3$dich_reneg[dataV3$V3_fav==1], dataV3$dich_reneg[dataV3$V3_bal==1])
fav_eq$p.value #p=0.21

# Change in percent supporting renegotiation from balanced to unfavorable
eq_unfav <- t.test(dataV3$dich_reneg[dataV3$V3_bal==1], dataV3$dich_reneg[dataV3$V3_unfav==1])
eq_unfav$estimate[1] - eq_unfav$estimate[2] #-6%
eq_unfav$p.value #p=0.01

# Analysis for Figure 6
##
## Means by each treatment combination for fairness of agreement
##

# Create subset of those who answered fairness dependent variable
dataV3 <- subset(data, is.na(data$V3_fair)==FALSE)

# Calculate means and confidence intervals based on equity (productivity) and trade balance treatment combinations

V3_fav_prod <- mean(dataV3$V3_fair[dataV3$V3_fav==1 & dataV3$V3_prod==1])
se.V3_fav_prod <- sqrt(var(dataV3$V3_fair[dataV3$V3_fav==1 & dataV3$V3_prod==1])/length(dataV3$V3_fair[dataV3$V3_fav==1 & dataV3$V3_prod==1])) 
ci.up.V3_fav_prod <- V3_fav_prod + se.V3_fav_prod * c.value   
ci.low.V3_fav_prod <- V3_fav_prod  - se.V3_fav_prod * c.value 

V3_fav_unprod <- mean(dataV3$V3_fair[dataV3$V3_fav==1 & dataV3$V3_unprod==1])
se.V3_fav_unprod <- sqrt(var(dataV3$V3_fair[dataV3$V3_fav==1 & dataV3$V3_unprod==1])/length(dataV3$V3_fair[dataV3$V3_fav==1 & dataV3$V3_unprod==1])) 
ci.up.V3_fav_unprod <- V3_fav_unprod + se.V3_fav_unprod * c.value   
ci.low.V3_fav_unprod <- V3_fav_unprod  - se.V3_fav_unprod * c.value 

V3_fav_equal <- mean(dataV3$V3_fair[dataV3$V3_fav==1 & dataV3$V3_equal==1])
se.V3_fav_equal <- sqrt(var(dataV3$V3_fair[dataV3$V3_fav==1 & dataV3$V3_equal==1])/length(dataV3$V3_fair[dataV3$V3_fav==1 & dataV3$V3_equal==1])) 
ci.up.V3_fav_equal <- V3_fav_equal + se.V3_fav_equal * c.value   
ci.low.V3_fav_equal <- V3_fav_equal  - se.V3_fav_equal * c.value 

V3_unfav_prod <- mean(dataV3$V3_fair[dataV3$V3_unfav==1 & dataV3$V3_prod==1])
se.V3_unfav_prod <- sqrt(var(dataV3$V3_fair[dataV3$V3_unfav==1 & dataV3$V3_prod==1])/length(dataV3$V3_fair[dataV3$V3_unfav==1 & dataV3$V3_prod==1])) 
ci.up.V3_unfav_prod <- V3_unfav_prod + se.V3_unfav_prod * c.value   
ci.low.V3_unfav_prod <- V3_unfav_prod  - se.V3_unfav_prod * c.value 

V3_unfav_unprod <- mean(dataV3$V3_fair[dataV3$V3_unfav==1 & dataV3$V3_unprod==1])
se.V3_unfav_unprod <- sqrt(var(dataV3$V3_fair[dataV3$V3_unfav==1 & dataV3$V3_unprod==1])/length(dataV3$V3_fair[dataV3$V3_unfav==1 & dataV3$V3_unprod==1])) 
ci.up.V3_unfav_unprod <- V3_unfav_unprod + se.V3_unfav_unprod * c.value   
ci.low.V3_unfav_unprod <- V3_unfav_unprod  - se.V3_unfav_unprod * c.value 

V3_unfav_equal <- mean(dataV3$V3_fair[dataV3$V3_unfav==1 & dataV3$V3_equal==1])
se.V3_unfav_equal <- sqrt(var(dataV3$V3_fair[dataV3$V3_unfav==1 & dataV3$V3_equal==1])/length(dataV3$V3_fair[dataV3$V3_unfav==1 & dataV3$V3_equal==1])) 
ci.up.V3_unfav_equal <- V3_unfav_equal + se.V3_unfav_equal * c.value   
ci.low.V3_unfav_equal <- V3_unfav_equal  - se.V3_unfav_equal * c.value 

V3_bal_prod <- mean(dataV3$V3_fair[dataV3$V3_bal==1 & dataV3$V3_prod==1])
se.V3_bal_prod <- sqrt(var(dataV3$V3_fair[dataV3$V3_bal==1 & dataV3$V3_prod==1])/length(dataV3$V3_fair[dataV3$V3_bal==1 & dataV3$V3_prod==1])) 
ci.up.V3_bal_prod <- V3_bal_prod + se.V3_bal_prod * c.value   
ci.low.V3_bal_prod <- V3_bal_prod  - se.V3_bal_prod * c.value 

V3_bal_unprod <- mean(dataV3$V3_fair[dataV3$V3_bal==1 & dataV3$V3_unprod==1])
se.V3_bal_unprod <- sqrt(var(dataV3$V3_fair[dataV3$V3_bal==1 & dataV3$V3_unprod==1])/length(dataV3$V3_fair[dataV3$V3_bal==1 & dataV3$V3_unprod==1])) 
ci.up.V3_bal_unprod <- V3_bal_unprod + se.V3_bal_unprod * c.value   
ci.low.V3_bal_unprod <- V3_bal_unprod  - se.V3_bal_unprod * c.value 

V3_bal_equal <- mean(dataV3$V3_fair[dataV3$V3_bal==1 & dataV3$V3_equal==1])
se.V3_bal_equal <- sqrt(var(dataV3$V3_fair[dataV3$V3_bal==1 & dataV3$V3_equal==1])/length(dataV3$V3_fair[dataV3$V3_bal==1 & dataV3$V3_equal==1])) 
ci.up.V3_bal_equal <- V3_bal_equal + se.V3_bal_equal * c.value   
ci.low.V3_bal_equal <- V3_bal_equal  - se.V3_bal_equal * c.value 

V3.means.f <- c(V3_bal_prod, V3_bal_unprod, V3_bal_equal, V3_unfav_prod, V3_unfav_unprod, V3_unfav_equal, V3_fav_prod, V3_fav_unprod, V3_fav_equal)
V3.ci.up.f <- c(ci.up.V3_bal_prod, ci.up.V3_bal_unprod, ci.up.V3_bal_equal, ci.up.V3_unfav_prod, ci.up.V3_unfav_unprod, ci.up.V3_unfav_equal, ci.up.V3_fav_prod, ci.up.V3_fav_unprod, ci.up.V3_fav_equal)
V3.ci.low.f <- c(ci.low.V3_bal_prod, ci.low.V3_bal_unprod, ci.low.V3_bal_equal, ci.low.V3_unfav_prod, ci.low.V3_unfav_unprod, ci.low.V3_unfav_equal, ci.low.V3_fav_prod, ci.low.V3_fav_unprod, ci.low.V3_fav_equal)
V3.names.f <- c("Bal. \n Prod.", "Bal. \n Unprod.", "Bal. \n Equal", "Unfav. \n Prod.", "Unfav. \n Unprod", "Unfav. \n Equal", "Fav. \n Prod.", "Fav. \n Unprod." , "Fav. \n Equal")

## Generate Figure 6 (adjust width/height as necessary)
par(mar=c(6,5,1,2))
plot(1:9, V3.means.f , main = "", 
     ylab = "", xlab = "", cex.lab=2,
     ylim = c(-.45, .55), xlim = c(0.75, 9.25), pch = 19, type = "p", xaxt="n", cex=1)
mtext("Fairness of Agreement", side=2, at=0.05, cex=1.4, line=2.5)
abline(h=0, col="blue", lty=2)
for(i in 1:length(V3.means.f)){
  lines(c(i,i), c(V3.ci.up.f[i], V3.ci.low.f[i]))
  mtext(V3.names.f[i], side=1, at=i, 3.12, cex=1.4)
}

# Analysis for quantities of interest cited in the tex in discussion of Figure 6

# Compare the equitable/favorable condition to the inequitable/favorable conditions
FavProd_FavUnprod <- t.test(dataV3$V3_fair[dataV3$V3_fav==1 & dataV3$V3_prod==1], dataV3$V3_fair[dataV3$V3_fav==1 & dataV3$V3_unprod==1])
FavProd_FavUnprod$estimate[1] - FavProd_FavUnprod$estimate[2] #0.17
FavProd_FavUnprod$p.value #p=0.01

FavProd_FavEqual <- t.test(dataV3$V3_fair[dataV3$V3_fav==1 & dataV3$V3_prod==1], dataV3$V3_fair[dataV3$V3_fav==1 & dataV3$V3_equal==1])
FavProd_FavEqual$estimate[1] - FavProd_FavEqual$estimate[2] #0.12
FavProd_FavEqual$p.value #p=0.09

# Compare the equitable/balanced condition to the inequitable/balanced conditions
BalEqual_BalProd <- t.test(dataV3$V3_fair[dataV3$V3_bal==1 & dataV3$V3_equal==1], dataV3$V3_fair[dataV3$V3_bal==1 & dataV3$V3_prod==1])
BalEqual_BalProd$estimate[1] - BalEqual_BalProd$estimate[2] #0.28
BalEqual_BalProd$p.value #p<0.01

BalEqual_BalUnprod <- t.test(dataV3$V3_fair[dataV3$V3_bal==1 & dataV3$V3_equal==1], dataV3$V3_fair[dataV3$V3_bal==1 & dataV3$V3_unprod==1])
BalEqual_BalUnprod$estimate[1] - BalEqual_BalUnprod$estimate[2] #0.30
BalEqual_BalUnprod$p.value #p<0.01


#
#
# The following analysis is for the supplementary appendix
#
#

# Generate quantities for Appendix Tables 1 and 2: Study Demographics
# create income brackets
data$Income_0_50 <- ifelse(data$income==1 | data$income==2, 1, 0)
data$Income_50_100 <- ifelse(data$income==3 | data$income==4, 1, 0)
data$Income_100_150 <- ifelse(data$income==5 | data$income==6, 1, 0)
data$Income_150up <- ifelse(data$income==7 | data$income==8 |data$income== 9, 1, 0)

# create party variables
data$Democrat <- ifelse(data$party==1, 1, 0)
data$Republican <- ifelse(data$party==2, 1, 0)
data$Independent <- ifelse(data$party==3, 1, 0)

#create education variable
data$CollegeDegree <- ifelse(data$edu>3, 1, 0)

#create age brackets
data$Age18_29 <- ifelse(data$age<30 & data$age>17, 1, 0 )
data$Age30_44 <- ifelse(data$age<45 & data$age>29, 1, 0 )
data$Age45_59 <- ifelse(data$age<60 & data$age>44, 1, 0 )
data$Age60up <- ifelse(data$age>59, 1, 0 )

#age variable
table(data$age)
data$age0 <- ifelse(data$age>= 18 & data$age <= 29,1,0)
data$age1 <- ifelse(data$age>= 30 & data$age <= 44,1,0)
data$age2 <- ifelse(data$age >= 45 & data$age <= 64,1,0)
data$age3 <- ifelse(data$age >=65,1,0)
data$agerange <- ifelse(data$age0==1, 1,ifelse(data$age1==1, 2,ifelse(data$age2==1, 3, ifelse(data$age3==1, 4, NA))))

# V1 - subset to those who answered support dependent variable
dataV1 <- subset(data, is.na(data$V1_supp)==FALSE )

# V3 - subset to those who answered support dependent variable
dataV3 <- subset(data, is.na(data$V3_App_reneg)==FALSE ) 


# Produce .tex code for Table 1
stargazer(dataV1[c("Age18_29", "Age30_44", "Age45_59", "Age60up", "women",  "Income_0_50", "Income_50_100", "Income_100_150", "Income_150up", "Democrat", "Republican", "Independent")], type = "latex", title="Study 1 Demographics", 
          summary.stat = c("mean"), digits = 2)

# Produce .tex code for Table 2
stargazer(dataV3[c("Age18_29", "Age30_44", "Age45_59", "Age60up", "women","Income_0_50", "Income_50_100", "Income_100_150", "Income_150up", "Democrat", "Republican", "Independent")], type = "latex", title="Study 2 Demographics", 
          summary.stat = c("mean"), digits = 2)

# Conduct balance tests for Appendix tables 3, 4, and 5

# Balance tests to check whether demographics are corellated with treatment assignment
Study1equal <- lm(V1_equal ~ age + edu2 + income  + women, data=dataV1)
Study1fav <- lm(V1_fav ~ age + edu2 + income + women, data=dataV1)
Study1unfav <- lm(V1_unfav ~ age + edu2 + income +  women, data=dataV1)

Study2bal <- lm(V3_bal ~ age + edu2 + income +  women, data=dataV3)
Study2fav <- lm(V3_fav ~ age + edu2 + income + women, data=dataV3)
Study2unfav <- lm(V3_unfav ~ age + edu2 + income +  women, data=dataV3)

Study2equal <- lm(V3_equal ~ age + edu2 + income + women, data=dataV3)
Study2prod <- lm(V3_prod ~ age + edu2 + income + women, data=dataV3)
Study2unprod <- lm(V3_unprod ~ age + edu2 + income +  women, data=dataV3)

# tex output for table 3
stargazer(Study1equal, Study1fav, Study1unfav,  keep.stat = c("n"), no.space=TRUE)
# tex output for table 4
stargazer(Study2bal, Study2fav, Study2unfav,  keep.stat = c("n"), no.space=TRUE )
# tex output for table 5
stargazer(Study2equal, Study2prod, Study2unprod,  keep.stat = c("n"), no.space=TRUE )


# Analysis for Appendix Figure 1, panels a, b, and c

# Mediation analysis comparing equal versus favorable agreements
# next lines before the "sensititivy test" repeats earlier code so you don't have to run all the preceeding code)
# subset to equal and favorable treatment
set.seed(43216)
dataV1eq.fav <- subset(dataV1[dataV1$V1_equal==1 | dataV1$V1_fav==1, ])

## fairness of agreement as a mediator
fair2.med <- lm(V1_fair ~ V1_equal + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1eq.fav)
fair2.dv <- lm(V1_supp ~ V1_fair + V1_equal + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1eq.fav)
fair2.med2 <- mediate(fair2.med, fair2.dv, treat="V1_equal", mediator="V1_fair", sims=1500)

## Sensitivity test for panel "a"
sens.fair.med <- medsens(fair.med2)
par.orig <- par(mfrow = c(2,2))
plot(sens.fair.med, main="", ylim=c(-1,1))

# Mediation analysis comparing equal versus unfavorable agreements 
# next lines before the "sensititivy test" repeats earlier code so you don't have to run all the preceeding code)
# subset to equal and unfavorable treatment
set.seed(43216)
dataV1eq.unfav <- subset(dataV1[dataV1$V1_equal==1 | dataV1$V1_unfav==1, ])

## fairness of agreement as a mediator
fair.med <- lm(V1_fair ~ V1_equal + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1eq.unfav)
fair.dv <- lm(V1_supp ~ V1_fair + V1_equal + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1eq.unfav)
fair.med2 <- mediate(fair.med, fair.dv, treat="V1_equal", mediator="V1_fair", sims=1500)

## Sensitivity test for panel "b"
sens.fair.med <- medsens(fair.med2)
par.orig <- par(mfrow = c(2,2))
plot(sens.fair.med, main="", ylim=c(-1,1))

# Mediation analysis comparing favorable versus unfavorable agreements
# next lines before the "sensititivy test" repeats earlier code so you don't have to run all the preceeding code)
# subset to favorable and unfavorable treatment
set.seed(43216)
dataV1fav.unfav <- subset(dataV1[dataV1$V1_unfav==1 | dataV1$V1_fav==1, ])

## fairness of agreement as a mediator
fair3.med <- lm(V1_fair ~ V1_fav + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1fav.unfav)
fair3.dv <- lm(V1_supp ~ V1_fair + V1_fav + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1fav.unfav)
fair3.med2 <- mediate(fair3.med, fair3.dv, treat="V1_fav", mediator="V1_fair", sims=1500)

## Sensitivity test for panel "c"
sens.fair.med <- medsens(fair.med2)
par.orig <- par(mfrow = c(2,2))
plot(sens.fair.med, main="", ylim=c(-1,1))

# Analysis for Appendix Figure 2

# Create an individual-level equality measure from the three equality questions
data$equal <-c()
for(i in 1:length(data$fair1_3)){
  data$equal[i] <- sum(data$fair1_3[i], data$fair2_3[i], data$fair3_3[i], na.rm = FALSE)
}

# Rescale for ease of interpretability
data$equal <- rescale(data$equal)


# Use quantiles to determine cut-points for individual-level measure low versus high equality for Trade Balance study
dataV1 <- subset(data, is.na(data$V1_supp)==FALSE ) 
quantile(dataV1$equal, c(0.25, 0.75), na.rm=TRUE) #25%=0.6 75%=0.8

#create low vs high equality indicator
dataV1$HighEqual <- ifelse(dataV1$equal > 0.79, 1, ifelse(dataV1$equal < 0.61, 0, NA))

# subset data into low and high equality
dataV1.LowEqual <- subset(dataV1[dataV1$equal < 0.61 & is.na(dataV1$equal)==FALSE, ])
dataV1.HighEqual <- subset(dataV1[dataV1$equal > 0.79 & is.na(dataV1$equal)==FALSE, ])
# subset to those who answered fairness question
dataV1fair.LowEqual <- subset(dataV1.LowEqual[is.na(dataV1.LowEqual$V1_fair)==FALSE, ])
dataV1fair.HighEqual <- subset(dataV1.HighEqual[is.na(dataV1.HighEqual$V1_fair)==FALSE, ])

# Generate fairness means and confidence intervals for treatmetns based on low/high equality
#Low Equality
V1_equal <- mean(dataV1fair.LowEqual$V1_fair[dataV1fair.LowEqual$V1_equal==1])
se.V1_equal <- sqrt(var(dataV1fair.LowEqual$V1_fair[dataV1fair.LowEqual$V1_equal==1])/length(dataV1fair.LowEqual$V1_fair[dataV1fair.LowEqual$V1_equal==1])) 
ci.up.V1_equal <- V1_equal + se.V1_equal * c.value   
ci.low.V1_equal <- V1_equal  - se.V1_equal * c.value 

V1_unfav <- mean(dataV1fair.LowEqual$V1_fair[dataV1fair.LowEqual$V1_unfav==1])
se.V1_unfav <- sqrt(var(dataV1fair.LowEqual$V1_fair[dataV1fair.LowEqual$V1_unfav==1])/length(dataV1fair.LowEqual$V1_fair[dataV1fair.LowEqual$V1_unfav==1])) 
ci.up.V1_unfav <- V1_unfav + se.V1_unfav * c.value   
ci.low.V1_unfav <- V1_unfav  - se.V1_unfav * c.value

V1_fav <- mean(dataV1fair.LowEqual$V1_fair[dataV1fair.LowEqual$V1_fav==1])
se.V1_fav <- sqrt(var(dataV1fair.LowEqual$V1_fair[dataV1fair.LowEqual$V1_fav==1])/length(dataV1fair.LowEqual$V1_fair[dataV1fair.LowEqual$V1_fav==1])) 
ci.up.V1_fav <- V1_fav + se.V1_fav * c.value   
ci.low.V1_fav <- V1_fav  - se.V1_fav * c.value

means.LowEqual <- c(V1_equal, V1_unfav, V1_fav)
low.ci.LowEqual <- c(ci.low.V1_equal, ci.low.V1_unfav, ci.low.V1_fav)
up.ci.LowEqual <- c(ci.up.V1_equal, ci.up.V1_unfav, ci.up.V1_fav)
V1.names <- c("Equal", "Unfavorable", "Favorable")

# High Equality
V1_equal <- mean(dataV1fair.HighEqual$V1_fair[dataV1fair.HighEqual$V1_equal==1])
se.V1_equal <- sqrt(var(dataV1fair.HighEqual$V1_fair[dataV1fair.HighEqual$V1_equal==1])/length(dataV1fair.HighEqual$V1_fair[dataV1fair.HighEqual$V1_equal==1])) 
ci.up.V1_equal <- V1_equal + se.V1_equal * c.value   
ci.low.V1_equal <- V1_equal  - se.V1_equal * c.value 

V1_unfav <- mean(dataV1fair.HighEqual$V1_fair[dataV1fair.HighEqual$V1_unfav==1])
se.V1_unfav <- sqrt(var(dataV1fair.HighEqual$V1_fair[dataV1fair.HighEqual$V1_unfav==1])/length(dataV1fair.HighEqual$V1_fair[dataV1fair.HighEqual$V1_unfav==1])) 
ci.up.V1_unfav <- V1_unfav + se.V1_unfav * c.value   
ci.low.V1_unfav <- V1_unfav  - se.V1_unfav * c.value

V1_fav <- mean(dataV1fair.HighEqual$V1_fair[dataV1fair.HighEqual$V1_fav==1])
se.V1_fav <- sqrt(var(dataV1fair.HighEqual$V1_fair[dataV1fair.HighEqual$V1_fav==1])/length(dataV1fair.HighEqual$V1_fair[dataV1fair.HighEqual$V1_fav==1])) 
ci.up.V1_fav <- V1_fav + se.V1_fav * c.value   
ci.low.V1_fav <- V1_fav  - se.V1_fav * c.value

means.HighEqual <- c(V1_equal, V1_unfav, V1_fav)
low.ci.HighEqual <- c(ci.low.V1_equal, ci.low.V1_unfav, ci.low.V1_fav)
up.ci.HighEqual <- c(ci.up.V1_equal, ci.up.V1_unfav, ci.up.V1_fav)

## Generate Appendix Figure 2 (adjust width/height as necessary)
par(mar=c(3,5,1,2))
plot(1:3, means.HighEqual , main = "", 
     ylab = "", xlab = "", cex.lab=2,
     ylim = c(-.45, .85), xlim = c(0.75, 3.25), pch = 19, type = "p", xaxt="n", cex=1.5)
mtext("Average Fairness Score", side=2, at=.2, cex=1.9, line=2.5)
abline(h=0, col="blue", lty=2)
for(i in 1:length(low.ci.HighEqual)){
  lines(c(i,i), c(up.ci.HighEqual[i], low.ci.HighEqual[i]))
  mtext(V1.names[i], side=1, at=i+.05, 1.12,cex=1.7)
}
points(1.15:3.15, means.LowEqual, pch=22, cex=1.5)
for(i in 1:length(means.LowEqual)){
  lines(c(i+.15,i+.15), c(up.ci.LowEqual[i], low.ci.LowEqual[i]))
}
text(2.97, .85, labels = "High Equality", cex=1.5)
text(2.966, .75, labels = "Low Equality", cex=1.5)
points(2.52, .85, pch=19, cex=1.5)
points(2.52, .75, pch=22, cex=1.5)

# Analysis for quantities of interest cited in the appendix text discussing Figure 2

## Effects of low/high equality in the equal treatment
high.low.bal <- t.test(dataV1$V1_fair[dataV1$HighEqual==1 & dataV1$V1_equal==1], dataV1$V1_fair[dataV1$HighEqual==0 & dataV1$V1_equal==1])
high.low.bal$estimate[1] - high.low.bal$estimate[2] #0.36
high.low.bal$p.value #p<0.01

# Test the treatment effect on low versus high equality respondents using an interaction effect 

# Differential treatment effect of equal vs unfavorable for high versus low equality respondents
Equality.equal.fav <- lm(V1_fair ~ V1_equal*HighEqual, data = dataV1[dataV1$V1_unfav != 1, ])
summary(Equality.equal.fav)# interaction 0.28, p<0.01
# Differential treatment effect of equal vs favorable for high versus low equality respondents
Equality.equal.unfav  <- lm(V1_fair ~ V1_equal*HighEqual, data = dataV1[dataV1$V1_fav != 1, ])
summary(Equality.equal.unfav)# interaction 0.22, p<0.03

# Analysis for appendix section 1.6 discussed in the appendix text

# Change in support for renegotiation from balanced/equal to balanced/productive treatments
BalEqual_BalProd <- t.test(dataV3$V3_App_reneg[dataV3$V3_bal==1 & dataV3$V3_equal==1], dataV3$V3_App_reneg[dataV3$V3_bal==1 & dataV3$V3_prod==1])
BalEqual_BalProd$estimate[1] - BalEqual_BalProd$estimate[2] #-0.01
BalEqual_BalProd$p.value #p=0.90
# Change in support for renegotiation from balanced/equal to balanced/unproductive treatments
BalEqual_BalUnprod <- t.test(dataV3$V3_App_reneg[dataV3$V3_bal==1 & dataV3$V3_equal==1], dataV3$V3_App_reneg[dataV3$V3_bal==1 & dataV3$V3_unprod==1])
BalEqual_BalUnprod$estimate[1] - BalEqual_BalUnprod$estimate[2] #0.03
BalEqual_BalUnprod$p.value #p=0.65



